% Simulate data from an SVAR(0) and characterise likelihoods under a 
% narrative restriction on the historical decomposition.

clear variables
close all
clc

rng(23032021); % Set random seed

T = 3; % Length of time series
n = 1000; % Number of gridpoints over which to evaluate likelihood
k = 1; % Period in which first structural shock is positive and main contributor
Keps = 1e6; % Monte Carlo sample size for estimating probabilities

% Set true values of structural parameters.
A0 = [1 0.5; 0.2 1.2];
% Matrix of impulse responses.
B = inv(A0);

% Compute reduced-form parameters
% Lower-triangular Cholesky factorisation of innovation covariance matrix.
Sigmatr = chol(inv(A0)*inv(A0)','lower'); 
Sigma = Sigmatr*Sigmatr'; % Covariance matrix of VAR innovations.
sig11 = Sigmatr(1,1);
sig21 = Sigmatr(2,1);
sig22 = Sigmatr(2,2);

% Construct grid for theta over which to evaluate likelihood.
thetaGrid = linspace(-pi+eps,pi-eps,n)';

%% Compute ex ante probability that restrictions are satisfied.
% Draw structural shocks from bivariate standard normal distribution.
epsk = randn([Keps,2]);

PrEps = zeros(n,1);

for ii = 1:n

    % Compute ex ante probability that structural shocks satisfy NR at
    % each value of theta on grid.
    PrEps(ii) = compPrEps(Sigmatr,thetaGrid(ii),epsk);

end

clear epsk

% Plot probability that restrictions are satisfied as a function of theta.
figure;
plot(thetaGrid,PrEps,'color',[0 0 204]./255,'LineWidth',2);
xlim([-pi pi]);
xticks([-pi -pi/2 0 pi/2 pi]);
xticklabels({'-\pi','-\pi/2','0','\pi/2','\pi'});
xlabel('\theta');
title('Probability Narrative Restrictions Satisfied');
Ax = gca;
Ax.FontSize = 14;
cd('Figures');
print('HistoricalDecompositionProbability','-depsc');
print('HistoricalDecompositionProbability','-dpng');
cd ..

%% Compute likelihoods under a particular realisation of the data.
% Draw a time series of structural shocks such that the first structural 
% shock in the kth period is positive and the most important contributor 
% to the observed change in the first variable.
epst = randn([T,2]);
flag = 0;

while flag == 0

    epst(k,1) = randn;
    flag = ((abs(B(1,1)*epst(k,1)) >= abs(B(1,2)*epst(k,2))) & ...
        epst(k,1) >= 0);

end

% Use shocks and structural parameters to generate data.
yt = (A0\epst')';

conLik = zeros(n,1);
unconLik = zeros(n,1);

for ii = 1:n % For each value of theta

    theta = thetaGrid(ii);
    
    % Value of first structural shock given value of theta and realisation
    % of data in period k.
    eps1k = (1/(sig11*sig22))*(sig22*yt(k,1)*cos(theta) ...
        + (sig11*yt(k,2)-sig21*yt(k,1))*sin(theta));
    
    % Check whether narrative restriction and sign normalisations are 
    % satisfied at value of theta.
    
    if cos(theta) >= 0 && sig22*cos(theta) - sig21*sin(theta) >= 0 ...
            && eps1k >= 0
        
        % Value of second structural shock given value of theta and 
        % realisation of data in period k.        
        eps2k = (1/(sig11*sig22))*((sig11*yt(k,2)-sig21*yt(k,1))*cos(theta) ...
            - sig22*yt(k,1)*sin(theta));
  
        if abs(sig11*cos(theta)*eps1k) >= abs(-sig11*sin(theta)*eps2k)
            
            % Compute unconditional likelihood.
            unconLik(ii) = exp(- T*log(2*pi) - (T/2)*log(det(Sigma)) ...
                - 0.5*trace((yt/Sigma)*yt'));            
            
            % Compute conditional likelihood.
            conLik(ii) = unconLik(ii)/PrEps(ii);
            
        end
        
    elseif cos(theta) < 0 && sig22*cos(theta) - sig21*sin(theta) >= 0 ...
            && eps1k >= 0
    
        % Value of second structural shock given value of theta and 
        % realisation of data in period k. 
        eps2k = (1/(sig11*sig22))*(sig22*yt(k,1)*sin(theta) ...
            -(sig21*yt(k,1) + sig11*yt(k,2))*cos(theta));
        
        if abs(sig11*cos(theta)*eps1k) >= abs(sig11*sin(theta)*eps2k)
            
            % Compute unconditional likelihood.
            unconLik(ii) = exp(- T*log(2*pi) - (T/2)*log(det(Sigma)) ...
                - 0.5*trace((yt/Sigma)*yt'));            
            
            % Compute conditional likelihood.
            conLik(ii) = unconLik(ii)/PrEps(ii);       

        end
        
    end
     
end

% Plot conditional and unconditional likelihoods.
figure;
plot(thetaGrid,conLik,'color',[0 0 204]./255,'LineWidth',2);
hold on;
plot(thetaGrid,unconLik,'color',[255 0 0]./255,'LineWidth',2,'LineStyle','--');
xlim([-pi pi]);
xticks([-pi -pi/2 0 pi/2 pi]);
xticklabels({'-\pi','-\pi/2','0','\pi/2','\pi'});
xlabel('\theta');
title('Likelihood');
Ax = gca;
Ax.FontSize = 14;
legend('Conditional','Unconditional','Location','NorthWest');
legend boxoff;
cd('Figures');
print('HistoricalDecompositionLikelihood','-depsc');
print('HistoricalDecompositionLikelihood','-dpng');
cd ..

cd('Results');
save('HistoricalDecompositionIllustration_results.mat');
cd ..

%%
function prob = compPrEps(Sigmatr,theta,epsk)
% Compute probability that shocks satisfy narrative restrictions.
% Setting probability to zero when sign normalisation on diagonal elements
% of A_0 is violated.
% Inputs:
% - Sigmatr: lower-triangular Cholesky factor of Sigma.
% - theta: structural parameter
% - epsk: matrix containing Monte Carlo sample of structural shocks

sig11 = Sigmatr(1,1);
sig21 = Sigmatr(2,1);
sig22 = Sigmatr(2,2);

if cos(theta) >= 0 && sig22*cos(theta) - sig21*sin(theta) >= 0
    
    prob = mean((abs(sig11*cos(theta)*epsk(:,1)) >= ...
        abs(-sig11*sin(theta)*epsk(:,2))) & epsk(:,1) >= 0);
    
elseif cos(theta) < 0 && sig22*cos(theta) - sig21*sin(theta) >= 0
    
    prob = mean((abs(sig11*cos(theta)*epsk(:,1)) >= ...
        abs(sig11*sin(theta)*epsk(:,2))) & epsk(:,1) >= 0);
    
else 
    
    prob = 0;
    
end

end